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■ Abstract. We consider the computation of averaged coefficients for the homogenization of 

' elliptic partial differential equations. In this problem, like in many multiscale problems, a large 

number of similar computations parametrized by the macroscopic scale is required at the microscopic 
scale. This is a framework very much adapted to model order reduction attempts. 

The purpose of this work is to show how the reduced-basis approach allows to speed up the 
computation of a large number of cell problems without any loss of precision. The essential com- 
' ponents of this reduced-basis approach are the a posteriori error estimation, which provides sharp 

error bounds for the outputs of interest, and an approximation process divided into offline and 
online stages, which decouples the generation of the approximation space and its use for Galerkin 
' projections. 
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1 Introduction 

In this work, we study the numerical homogenization of linear scalar elliptic partial 
differential equations (PDEs) such as those encountered in the problems of thermal 
diffusion and electrical conduction. Oscillating test functions, also termed correctors, 
are computed through a reduced-basis (RB) approach for parametrized cell problems 
supplied with periodic boundary conditions. Numerical results have been obtained 
with some prototypical parametrizations of the oscillating coefficients and are shown 
in a two-dimensional case with one single varying rectangular inclusion inside rect- 
angular cells. The method applies to all numerical homogenization strategies that 
require to solve a large number of parametrized cell problems. 

In periodic homogenization, only one cell problem has to be solved in order to 
completely determine the homogenized coefficient (s) to be used in the homogenized 
(macroscopic) equation. In sharp contrast, non-periodic homogenization requires the 
solution of several cell problems (in fact, theoretically, an infinite number of them, 
and in practice, a large number). A homogenized coefficient is then approximated 
by some average over a large number of microscopic cells. Consequently, as opposed 
to the periodic case where the computation is light and exact, the non-periodic case 
asks for a computationally demanding and approximate-in-nature task. This is why 
the design of a fast and accurate numerical homogenization method is considered as 
an important issue for the treatment of non-periodic heterogeneous structures. The 
RB approach seems very well adapted to this framework. 

The article is organized as follows. In section 2, we give a detailed presentation 
of the setting of the problem. For the sake of consistency and the convenience of the 
reader, we also briefly outline the main relevant issues in homogenization and RB 
theories. In section 3, the RB approach for a parametrized cell problem is introduced 
and we notably derive a posteriori error bounds related to the convergence of the 
RB method in the homogenization context. Numerical results for the prototypical 
example of rectangular cells with one single rectanguler inclusion are presented in 
section 4. Possible extensions of our work are discussed in the final section. 

2 Setting of the problem, elements of homogenization 
theory and RB approach 

2.1 Formulation of the problem 

The mathematical problem under consideration throughout this article reads as 
follows. We are interested in the behaviour of a sequence of scalar functions that 
satisfy 



in a bounded open set Q C M", for a sequence of scalars e > 0. Of interest is the 
asymptotic limit of the sequence when e — > 0, along with approximations for n*^ 




(1) 
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when e is small. 

For the sake of simplicity, the scalar source term / is chosen in L'^{Q) and we 
supply equation ([1]) with the following boundary conditions on the smooth (say 
C°°-Lipschitz) boundary dQ = TdIJ Tjy of Q, 

But as a matter of fact, it is well known that the homogenization results are local 
in nature and do not depend on the boundary conditions, except for what regards 
error estimations close to the boundary. Nor do the homogenization results depend 
on the source term /. Hence the generality of the assumptions (BC) and / G L'^{Q), 
chosen here to give a precise mathematical frame to the numerical experiments. 

To fix ideas, the unknown could be thought of, either as a temperature or 
as an electric field in a macroscopic domain Q. The tensorial coefficients for A''{x) 
would respectively be thought of, either as temperature diffusivities or as electric 
conductivities. 

Next, let us define, for any e > 0, the family A"^ £ L°°{Q,^AaA,'rA) functions 
from Q to the set A^a^,7^ of uniformly positive definite nxn matrices (second order 
tensors) with uniformly positive definite inverses, that is, matrices satisfying, for 
all X S 0, 

<aA|n|2< A'{x)u-u,yueR"- (3) 
<7a|uP< I'(x)~^n-ti,Vn gR" . (4) 

Under such conditions, equations (HI)-© are well posed in the sense of Hadamard. 
For every e > 0, there exists a unique solution n*^ in H^^ (Q) = G H^{U),u |r£,= O} 
that continuously depends on /, 

\\u^HHn)<C{mf\\LHn) , (5) 

with some constant C(il) that is only function of fi. Moreover, the sequence of 
solutions n*^ is bounded in H^^{Q,), so that some subsequence e' weakly converges 
to a limit u* G i?p^(r2) when e' — > 0. We are specifically interested in estimating 
the behaviour of this weakly-convergent subsequence. 

In a typical frame for the homogenization theory, the coefficients A^ are as- 
sumed to oscillate very rapidly on account of numerous small heterogeneities in the 
domain For example, e typically denotes the ratio of the mean period for mi- 
croscopic fast oscillations of A^ divided by the mean period for macroscopic slow 
oscillations of A^ in fi. Moreover, it is usually assumed that macroscopic (macro) 
and microscopic (micro) scales "separate" when e is sufficiently small, which allows 
for the oscillating coefficients to be locally homogenized in the limit e ^ 0. 



2.2 General context for homogenization 

As announced above, this section [2^2] includes some basics of homogenization theory 
for linear scalar elliptic PDEs. The purpose of this summary is only to collect some 
elementary results for convenience. Readers familiar with the homogenization theory 
may then like to skip this section and proceed to section 12. 3[ which introduces the 
RB theory. 
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2.2.1 Abstract homogenization results 

The following abstract homogenization result is the basis for many studies that aim 
at computing a numerical approximation for when e is small [131 IT71 [21 ITTI [9] . 
It shows that, in the limit e — > 0, the small oscillating scale "disappears" from the 
macroscopic point of view ; that is, the microscopic and macroscopic behaviours 
asymptotically "separate". This implies that the limit problem is easier to solve 
than equation ([T]) for some small e, since the former does not require to resolve 
microscopic details. Moreover, a tractable approximation of when e is small 
enough can be computed from the asymptotic limit when e — > 0. 

More precisely, u* can be obtained as the solution to the H-limit equation for ([1]) 
(see equation ([8]) below). It is then an L^-approximation for u*" when e is small, as the 
asymptotic L^-limit of n*^ when e — > 0. Moreover, an improved if ^-approximation 
for when e is small can also be computed with u* after "correction" of the gradient 
Vn*. 

The homogenization of the sequence of equations ([T]) is the mathematical process 
which allows to define the H-limit equation and the approximation for u''. It is 
performed using the following abstract objects 



a sequence of n oscillating test functions £ H^(Q) such that, for every direction 
(ej)i<j<n of the ambient physical space R", we have zf^Xi in H^{Q) and 

-div(i^Vzl) = -div(i*ei) in H-\n) , 

a homogenized tensor A* defined by 

I^V4 ^i*ei in , (6) 

a subsequence u'^' of solutions for ([T]) that satisfies 



u* in H^n) 



(7) 



[ A' Vu' A^'Vu* in [L2(17)]" 
where n* is solution for the H-limit or homogenized equation 

- div{A*{x)Vu*{x)) = f{x), Vx G , (8) 
supplied with the boundary conditions (BC), 
and an asymptotic approximation for a subsequence e' of e that satisfies 

-u- (9) 

L2(n) 



1=1 



^—^ , (10) 

where diU* are the components of Vu* in each direction Cj. 

Note that the latter convergence result (Ho]) for Vu^' also holds in [Lj;^(17)]" 
if u* £ W^''^{i^). So, if u* £ ii^(r2), the corrector result states that u'' can be 
approximated with the following formula, 



n 

u 

i=l 



+ - Xi)diu'' + r, , (11) 
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where the remainder term converges strongly to zero in W^^^ (O). 

In a nutshell, the homogenization of the sequence of equations ([1]) has allowed 
to derive an abstract homogenized problem, ©-([H]), the solution u* of which can be 
corrected with (jlOp into an approximation of n*^ in the limit e — > 0. 

But we lack an explicit expression for the homogenized tensor A* to get an 
explicit asymptotic limit u*. That is why, though it is not required by the previous 
abstract theory, the scale separation in the behaviour of the oscillating coefficients 

is often assumed to be explicitly encoded, using some specific postulated form for 
A^. This allows to derive an explicit expression of the homogenized problem, and 
even an error estimate in terms of e for the correction error r^ in (llip . which allows 
to quantify the homogenization approximation error. 

2.2.2 The explicit two-scale homogenization 

To get explicit expressions for the homogenized problem, some particular depen- 
dence of the family A^ on the space variable x is often assumed, like in two-scale 
homogenization for instance. Namely, on account of the scale separation assumption 
and the local dependence of the homogenization process, one of the most common 
assumption is the local periodicity for A"^, which can be made precise as follows. 

It is assumed that tensors A^ are traces of functions of two coupled variables on 
the set locally defined by a fast microscopic variable e~^x linearly coupled with the 
slow macroscopic variable x in 

A'{x) = a(x,^'^ , (12) 

where, for any x £ i}, the function A{x, •) 

A{x,-) :yeM"^i(x,y) € M"^'^ 

is 1-periodic in each of the n directions (ei)i<j<n, which makes the local oscilla- 
tions completely determined when e — > 0. The domain Y = [0, 1]" of the pe- 
riodic pattern is called the cell and is identified with the n-dimensional torus. 
A(x, •) is said to be Y-periodic. Note that the properties of the tensors A^ imply 

AeL^{n,L°-{Y,Ma^,^^)). 

Now, under the assumption (jl2p of local periodicity, one possible manner to get 
explicit expressions for the homogenized problem is to perform a formal two-scale 
analysis with the following Ansatz 




where, for any x £ ^1, the functions Ui{x,-) are Y-periodic. The first two terms of 
the Ansatz ()13p are shown to coincide with the approximation (jlip for it^ [31[T]. 

Inserting the Ansatz ()13p into equation ([T|) gives the following explicit expressions 
for the objects previously defined by the abstract homogenization result: 

• the function uq = u*{x) does not depend on the fast variable e~^x and is the 

approximation for given by the convergence result Q, 

• the gradient VyUi{x, •) linearly depends on VxU*{x), 

ui (^x, -) = ^ diU*{x)wi i^x, + Mi(x), 
where (u)i(x, •))i<i<„ are n Y-periodic cell functions, 
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• the n cell functions Wi{x,-), parametrized by their macroscopic position x E 

are solutions to the following n cell problems, 

-diYy{A{x,y)-[ei + VyWi{x,y)])=0,yy£W , (14) 
and the correctors now read z| = Xj + ewi{x,x/e), 

• the entries ( A*{x)ij I of the homogenized matrix A* can be explicitly com- 

V ' / l<i,j<n 

puted with the cell functions Wi(x, •), 

^*ix)i,j = J^A{x,y)[ei + VyWi{x,y)]- Cj dy, (15) 

• the approximation for u"" is now tractable and writes 

n 

= u* + WidiU* + r^ , (16) 

i=l 

where, provided u* S W'^'°°{Q), the correction error can be estimated to 
locally scale as e (far enough from the boundary layer), and to globally scale 
as ^/€, 

W^eWm < Cie\\u*\\w2,oo(^^),\/uj (sQ, (17) 

\\re\\H^^(n) < C2\^\\u*\\w2,c.^n), (18) 

with constants Ci depending only on 0,. 

To sum up, the local periodicity assumption (fT2]) allows to completely determine 
the homogenized problem through explicit two-scale expressions. The derivation of 
the homogenized equation in the case of locally periodic coefficients serves as a basis 
for many numerical homogenization strategies. 

2.2.3 Numerical homogenization strategies 

Under local periodicity assumption ()12p . a two-scale explicit homogenization strat- 
egy for a sequence of linear scalar elliptic PDEs like ([T]) reads as follows in the frame 
of Finite-Element approximations for the scalar elliptic problems ()14p and (jH]) . 

Algorithm 1 (Two-scale homogenization strategy) To homogenize the sequence 
of PDEs (CP.- 

1. solve the parametrized cell problems at each point x £ where the value 
of A*{x) is necessary to compute the FE matrix of the homogenized problem 

m, 

2. store the functions Wi for the future computation of the approximation of 

3. assemble the FE matrix associated with the homogenized operator — div(A*V-), 
4- solve the macroscopic homogenized problem (0], 

5. build the approximation lllb]) for u'' with u* and Wi. 
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On the one hand, in many practical situations, it is very common to assume 
that the tensors satisfy assumption (jl2p . Indeed, in practice, is often known 
for some given e = cq only. So, the asymptotic structure A'' of the problem with 
oscillating coefficients in Q has to be constructed from the only member A^'^ . Now, 
to take advantage of the exact explicit expressions given by the two-scale analysis, it 
is preferable to build a family A^ that satisfies assumption (jl2p . when possible. As- 
sumption ()12|) then seems fully justified for many applications from the practitioner's 
point of view. Then, the main numerical difficulty of the two-scale homogenization 
is the first step, that is the accurate computation of a large number of cell functions. 
This is the main issue addressed in this article. 

On the other hand, for some applications where the heterogeneities are highly 
non-periodic, one may want to build the sequence A^ differently, or even skip the ex- 
plicit construction of the sequence A"". For example, the actual construction process 
of heterogeneities may suggest another sequence A'' for which the error estimation of 
the approximation would then be more precise and meaningful. Or it may seem 
too difficult to explicitly build such a sequence A^ that satisfies ()12p from the knowl- 
edge of some only. In such cases, many numerical homogenization strategies 
have been developped to treat the numerical homogenization of oscillating coeffi- 
cients that are not locally periodic. 

To our knowledge, most of the existing numerical homogenization strategies may 
be classified in one of the two following categories. They either rely on different space 
assumptions than local periodicity for the oscillating coefficients {e.g. reiterated 
homogenization [15], stochastic homogenization [6j, deformed periodic coefficients 
[8], stochastically deformed periodic coefficients [5]), and still allow to derive exact 
(but not always fully explicit) expressions for the homogenized equation and the 
error estimate of the approximation. 

Or, the numerical homogenization strategies are much coarser and only rely on 
the assumption that explicit scale separation allow for the behaviour of the oscillating 
coefficients at some small e to be numerically homogenized. Those strategies are 
then approximate-in-nature. They manage to approximate quite a large class of 
heterogeneous problems, but may be computationally very demanding. They may 
also lack sharp error estimates. Example are the Multiscale finite-element method 
(MsFEM) [ISIE], the Heterogeneous multiscale method (HMM) [9], or the recent 
variational approach for non- linear monotone elliptic operators proposed in [llj... 

Now, in any of the two previously described situations where the numerical ho- 
mogenization strategies require the computation of a large number of parametrized 
cell problems, the RB approach proposed here-after is likely to bring some addi- 
tional computational efficiency. As a matter of fact, most numerical approximate 
homogenization strategies are only slight modifications of the exact two-scale ho- 
mogenization strategy proposed above in the frame of local periodicity assumption 
(I12p . and they do require the computation of a large number of parametrized cell 
functions. For many mechanical applications, this owes to the assumed existence of 
a Representative volume elements (RVE), which leads to general cell problems at 
each point x of the macro domain [10] . One simple example of a possibly approx- 
imate numerical homogenization strategy that requires the computation of a large 
number of parametrized cell functions is based on the following theorem, proved by 
Jikov et al. in ^4|. 

Theorem 1 Let A^ be a sequence of matrices in L°°(r2, TWq,^^) that defines a se- 
quence of linear scalar elliptic problems like |2|j. The H-limit of A'' is the homoge- 
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nized tensor A* . 

For any x G and e > 0, and any sufficiently small h > 0, let us define a 
sequence of locally periodic matrices (in the sense of (fT2]) ). 

AUx) = A'{x + h[e-^x]) , (19) 

where [e~^x] denotes the integer part ofe'^x. 

Then, for every 1 < i < n, there exists a unique sequence of periodic solutions 
w^'^{x, •) in the quotiented Soholev space H}^{Y)/M. of Y -periodic functions in H^{Y) 
that satisfy the n cell problems 

- div(i^(x + hy) ■ [ci + Vywl^\y)]) = 0, Vy G F (20) 

in the n-torus Y = [0, 1]". 

For each point x G and e > 0, we define a matrix A*^ made of elements 

■ = / ^'(^ + ■ + ^y^i'^i^^ y)] ■ [f^j + ^yW^-^ix, y)] dy 

for any {i,j) in {1,2, . . . Then, there exists a subsequence h' ^ such that 

lim limA*.,(x) = A*(x) . 

This theorem shows that, for any family A^, it is always possible to approximate 
the exact homogenized problem with the same explicit expressions than those ob- 
tained under the local periodicity assumption (I12p . after local periodization of A"^ 
hke in (fT9l) . 

So, considering the landscape for the homogenization theory as described above, 
among alternatives way of improving the numerical homogenization strategies, we 
choose here to concentrate on speeding up the numerical treatment of a large number 
of parametrized cell functions, rather than, for example, refining the approximations 
leading to explicit cell problems for larger class of oscillating coefficients. 

In the sequel, for the sake of simplicity, we assume that the sequence of tensors 
A^ satisfies assumption ()12p and apply the RB approach to the two-scale numerical 
homogenization strategy (Algorithm [T]) . Yet, the RB approach may apply as well 
with any numerical homogenization strategy that consists of first approximating 
A"^ by some sequence of tensors Aj^ like in (I19p . the latter leading to an explicit 
approximation for the homogenized problem after solving parametrized generalized 
cell problems like (|2Up . Let us now concentrate on decreasing as much as possible 
the computational cost of solving ()14p for many parameter values x G ^l. 

2.3 The reduced-basis method 

Two critical observations allow to think that an output-oriented model order re- 
duction technique like the RB method is likely to improve the repeated numerical 
treatment of parametrized cell problems ()14p . First, only outputs of the cell func- 
tions are required to solve the homogenized problem, and an a posteriori estimation 
gives sharp error bounds for those outputs. 

Second, as extensively discussed above, the numerous parametrized cell prob- 
lems arising from numerical homogenization strategies can be solved independently 
for each value of the parameter. Thus, a computational procedure based on an of- 
fline/online approach should naturally allow for a reduction of the computation time 
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in the limit of many queries. In particular, a large number of (and theoretically, an 
infinity of) parametrized cell problems occurs in the limit e — > of the homogeniza- 
tion strategies, in order to compute the homogenized problem ([8]) with non-periodic 
coefficients. And the number of homogenized problems to compute and solve can 
also be very large in practice, for instance in the frame of parameter estimation and 
optimization problems. 

These two observations motivate an RB approach for the parametrized cell prob- 
lems (I14p . which should significantly decrease the expense of computations in terms 
of CPU time for the homogenization problems where the offline stage is short com- 
pared to the online stage, or where the offline stage is not even an issue (like in 
real-time engineering problems for instance) [T^]. We are now going to introduce 
the basics of the reduced-basis method, well known to experts, who may want to 
directly proceed to the Section [3l 



2.3.1 The parametrized cell problem 

Let X be the quotiented space H^(Y)/'R of Y-periodic functions that belong to the 
Sobolev space H^{Y). The Hilbert space X is imbued with the H^{Y)-norm 

_l/2 

I^IIa" = ( / ■ 



induced by the inner product [u, v) = fy " v) ^ X x X- In the 

dual space X' of X, the dual norm is defined for any g £ X' hy 

II II Siv) 

WaWx' = sup ■ 



f6A' \\V\\x 

For any x £ ft, we define: 
a continuous and coercive bilinear form in X x- X parametrized hy x £ il,, 



a{u,v;x) = J A{x,y)Vu{y) ■Vv{y)dy, \f{u,v) e X X X , 

for which ua and are respectively coercivity and continuity constants, 
• and n continuous linear forms in X also parametrized by x E ^2, 

fi{v; x) = - j A{x, y)ei ■ Vv{y)dy, V?; G A" , 1 < i < n. 

Now, for any integer i, 1 < i < n, the i-th cell problem (|14|) for the cell functions 
Wi{x, •) rewrites in the following weak form: Find Wi{x, ■) £ X solution for 

a{wi{x,-),v;x) = fi{v]x) ,yv £ X , (21) 

where x £ il. plays the role of a parameter. 

We set Aii = {'Wi{x,y),x £ 0,} the solution subspace of the i-th cell problem 
(|2ip induced by the variations of x in 17, and 

TV 

M = {wi{x, y),x £n,l <i< N} = [jMi 

1=1 

the global solution subspace, that is the reunion of all solution subspaces for all cell 
problems. 
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Remark 1 Note at this point that Mi and Ai should be seen as spaces induced 
by the family of coefficients {A*{x,-)] , and not x. It is indeed always possible 

(and often useful) to use other explicit quantities than x as parameters to map Mi 
and M, provided that the variations of the parameters inside a given range of values 
induce the same family of coefficients and the same corresponding cell functions than 
X €Q. 

Besides, for the sake of simplicity in the presentation of the RB method, the 
tensor A{x, •) will be assumed symmetric in the following. Hence, in the computation 
of the homogenized tensor A*(x), the only interesting output for the n solutions 
Wi{x,-) is given by a symmetric matrix s of size n x n (somewhat similar to a 
compliance in the terminology of mechanics). The entries {sij)i<ij<n of the matrix 
s are given by 

Sijix) = -fj{wi{x,-);x) = j^A{x,y)Vwi{x,y) ■ ejdy . (22) 

But note that the RB approach still applies with non-symmetric tensors A{x,-) 
modulo slight modification^. 

The purpose of the RB method is to speed up the computation of a large number 
of solutions Wi{x,-) € X of (j2ip for many parameter values x G $7 while controlling 
the approximation error for the output s. 



2.3.2 Principle of the reduced-basis method 

The purpose of most order reduction techniques like the RB method is to solve 
the weak form of a PDE like (I2ip through a Galerkin projection method with a 
Hilbertian basis that is "adapted" to the solution subspace M. 

For instance, a Hilbertian basis that is adapted to the equations ()2ip when x G 
and 1 < z < n is an orthonormal family (^j)jgN (orthonormal with respect to the 
ambient inner product (•, •)a') such that 

• the ambient solution space X C spanj^j, j G N} is separable, 

• and, for a finite A^-dimensional subspace Xn = span{(^j, 1 < j < A*"} of X-, the 
Galerkin approximations Wiiy{x, •) G Xn for Wi{x, •) that satisfy, for any x in 
i7 and 1 < i < N, 

a{wiN{x,-),v;x) = fi{v;x) ,Mv G Xn , (23) 
are "sufficiently" close to Wi{x, •) for a given tolerable precision. 

But the previous definition is only vaguely stated until the "tolerable precision" is 
mathematically defined. 

One possible way of defining sufficient precision is to control the approxima- 
tion error for Wi{x,-) G X with the natural norm || • of the cell problem. The 
reduced-basis method rather proposes to control the approximation error for some 

^When the tensor A{x,-) is not symmetric, a dual problem, adjoint to the problem (|2ip . is 
introduced. The dual problem can be solved similarly to the "primal" problem (I21|l with an RB 
method, in a dual RB projection space. Last, the output should be rewritten like s plus an additional 
term that accounts for the residual error due to the RB projection of the equation (|2ip . This primal- 
dual approach is more extensively described in [18] for instance. 
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linear output like s, which is not very different in the present case where A{x,-) is 
symmetric as it will be made clearer in section I3.1i 

The RB method is based on the computation of a basis for a Galerkin projection 
space in X that is "adapted" to the cell problem (I14p in the sense of the minimiza- 
tion of the output approximation error. The approximation errors are made explicit 
through rigorous a posteriori estimates, which allow to a posteriori certify the ef- 
ficiency of the model order reduction, that is, the convergence of the RB method 
when the size of the Galerkin projection space increases. 

2.3.3 Practice of the reduced-basis method 

The RB method computes a basis for Xn from an approximation Mp of M, 

Mp = {wi{x, ■),Xk e V,! < i < n} , 

induced by a discrete sample 2? = {x^, l<k<p},pxn>N,oi parameter values 
distributed over the parameter space 0, (P C 0,). This is termed as the offline 
stage. Such a model order reduction is efficient if the tolerable precision for the 
output approximation error is reached with a A^-dimensional "adapted basis" when 
is mcuh smaller than M, where M is the number of degrees of freedom necessary 
for a generic numerical method, like the FE method, to reach the same precision. 
We call reduced basis such an "adapted basis" {S,j)i<j<nxN- 

Now, the RB treatment of ()2ip begins with the computation of a sample of cell 
functions that induces Aip. The cell functions of the finite space Mp should then 
be approximated before the model order reduction is possible. An accurate and 
generic numerical methoc|l, with a large number Af of degrees of freedom, like an 
FE method with a fine mesh for Y for instance, should be used at the beginning of 
the offline stage to compute an approximation for Aip. 

Then, the sample of solutions Aip from which the basis {(,j)i<j<nxN is built 
should be as representative as possible of A^. In the absence of information, 2? 
should be chosen arbitrarily. So, it is not only quite often impossible to choose 
a priori the "right" parameter sample D for Xn in order to minimize the output 
approximation error over every A^-dimensional vector subspace of M, but besides, 
V should not be too large so that the reduced-basis construction is fast compared 
to the online stage. Hence the necessity for a reliable a posteriori control of the RB 
approximation method, which allows to build fast a reduced basis {S,j)i<j<nxN, as 
it will be seen in section 13.21 

After building a reduced basis for the vector field Xn, the Galerkin projection 
method is applied to the weak form for (j2ip at any x £ ^l. That is, in this online 
stage, the previous reduced basis is assumed to span a vector field Xn sufficiently 
close to the solution manifold M of the parametrized PDF so that we can compute 
fast a sufficiently accurate Galerkin approximation in Xn for the solution of the 
parametrized PDF at any parameter value x £ i}. 

^In general non-symmetric cases, the approximation error for linear outputs like s can be ex- 
pressed as a product of two approximation errors, one for the parametrized solutions 'Wi{x, •) and 
another one for some dual quantity that is solution for the problem dual to (|2ip . But here, because 
of the symmetry and of the specific nature of the output, the so-called compliance in reference to 
mechanics, the dual problem is unnecessary and the approximation error for the output can be 
directly expressed as the square of the approximation error for the cell function. 

^Note that the time needed to compute a (possibly large) sample of p accurate FE approxi- 
mations can also be an issue, that can be dealt with by a pre-processing stage according to the 
parametrization, as it will be made clearer in section [32] 
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3 Reduced-basis approach for the cell problem 

The RB approach for equation (j2ip needs to a posteriori estimate the approximation 
error for Galerkin solutions of the cell problem. In a second stage, this allows for an 
a posteriori estimation of the approximation error on the output s. 

3.1 Error bounds for the cell problem 

The purpose of this section is to derive the error bound (|28p for the cell funtioncs, 
solutions of equation ()2ip . This allows to a posteriori estimate the approximation 
error for Galerkin solutions of the cell problem, and their outputs through the error 
bound (p3]) . To this end, let us introduce the linear operator : X ^ X so that, 
for any u £ X and x £ fi, 

{T^u,v)x = a{u,v;x), \/v e X ■ 

The existence of such an operator directly leans on the Riesz-Frechet representation 
Theorem in the Hilbert space X- 

For any 1 < i < N and x £ Q, the Galerkin approximation error 

\\wi{x,-) -WiN{x,-)\\x (24) 

can be bounded starting from the following equality, 

a{wi{x,-) - WiN{x,-),v;x) = fi{v;x) - a{wiN{x,y),v]x),Mv G X, (25) 

which is easily obtained by substraction of ([23]) from (f2T|) . 

Let us define the parametrized bilinear residual forms gi in X x X such that, for 
all parameter values x G and 1 < i < n, 

gi{u,v;x) = a{u,v;x) - fi{v;x),\/{u,v) e X x X ■ 

Then, equation ([25]) with v = Wi{x, •) — wnyi^x, •) allows to immediately derive the 
following estimates through the dual norm of the residual linear form for Wi{x,-) 
defined in X 

V gi{wiNix,-),v;x) . 
First, owing to the coercivity of the bilinear form a, we obtain the lower bound: 

aA\\wi{x,-) - WiN{x,-)\\x < \\gi{wiN{x,-),v;x)\\x' , (26) 

for the Galerkin approximation error. 

Second, in view of the continuity of the bilinear form a, we obtain the superior 
bound: 

\\gi{wiN{x,-),v;x)\\x' <7A^\\'^ii^r) -WiN{x,-)\\x , (27) 

for the Galerkin approximation error. 

Finally, note that it is possible to compute the dual norm of the linear form 

V ^ gi{wiN{x,-),v;x) = -a{wi{x,-) - WiNix,-),v]x) 

using the Riesz-Frechet representant T^{wi{x, •) — Wij^ix, •)) in the Hilbert space X, 
and that one can obtain numerical approximations for a a and 7^^, either using 
the spectral properties of the matrices resulting from the Galerkin projection in 
large generic solution spaces during the offline stage, or using properties of the 
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parametrization like in section 13.51 So, the Galerkin approximation error ()24p can 
be a posteriori bounded using estimations (j26|) and ([27|) . 

For X G n and 1 < z < A^, we define a posteriori estimators A]\f{wi{x, •)) for tlie 
Galerkin approximation errors (|24p . using the previous superior bounds, by 

^, / / ^^ ||a(w»(x, •) - W^Af (x, •), •; x) [1;^./ 

AAr(t(;j(x, •)) = • (28) 

aA 

The effectivities r]]\f{wi{x, •)) corresponding to the estimators Aj\f{wi{x, •)), 

r/AT U-i X, • = T, — 7 ^ 7]]—, (29) 

\\Wi{x,-) -WiN[X,-)\\x 

satisfy the following inequalities independently of N, 

i<miMx,-)) <—, (30) 

aA 

which shows the stability of the error estimator A]\f{wi{x, •)). 

Last, Galerkin approximations for the homogenized and the output matrix write 



^*Ni^)i,j = J A{x,y)[ei + VyWiN{x,y)] ■ Cj dy, (31) 
sfjix) = / A{x,y)^WiNix,y) ■ Cjdy . (32) 



JY 

The a posteriori superior bound A]\[{wi{x, •)) for the Galerkin approximation error 
()24p will now allow us to derive a simple superior bound for output approximation 
errors. Indeed, we have for any 1 < i, j < N and x £ Q, 

\ Sij{x) - sfj{x) \ = \ fj{Wi{x,-) - WiN{x,-);x) \ 

= I a{wj{x,-),Wi{x,-) - WiN{x,-)]x) I 

= I a{wj{x, •) - WjNix, ■),Wi{x, •) - WiN{x, ■);x) \ 

< aAANiwi{x,-))ANiwj{x,-)) 

since 'Wj{x, •) and Wi{x, •) are solutions for ([2T]) . and WiN{x, ■) is solution for ([23]) . 

We are finally in possession of an a posteriori superior bound A|^-^(x) for 
Galerkin approximations of the output Sij{x), 

AS t \ _ MMx,-) - 'WiNix,-),-,x)\\x'\\ai'Wjix,-) - WjNix,-),-;x)\\x' 

Numerical approximations for A|^. j^{x) will allow us to build fast a reduced basis for 
cell problems (f2T]) . Note that A*^- jsf{x) scales as the product A]\f{wi{x, ■))AN{wj{x, •)), 
hence the interest of model order reduction techniques for solutions Wi{x,-) without 
much loss of precision for output s{x). 

Remark 2 Note that for the output error bounds to scale like the square of the error 
bound for the cell functions, it has been essential to have the following orthogonality 
property for any x £ il., 

a{wi{x, ■) - WiN{x, WjNix, ■);x) =0 , VI < i,j < N. 

That is why we have chosen to build only one RB projection space Xn, spanned by 
all the parametrized cell functions Wi{xk,y) when 1 < i < n and Xk G D- Yet, note 
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that without this choice, the same scaling can still be obtained with n distinct RB 
projection spaces {XiN)i<i<n for each of the n solution subspaces Aii, provided one 
slightly modifies the definition of the output. Namely, another output matrix a and 
its RB approximation should then be defined, starting from s and s^ , by adding 
a residual error. Their entries read, for I < i, j < n, 

aij{x) = -fj{wi{x,-);x) + gi{wi{x, ■),Wj{x, ■); x) (34) 
crfj{x) = -fj{wiNix,-);x) + gi{wiN{x,-),WjNix,-);x) , (35) 

where a = s because the tensor A{x, •) is symmetric, and = s^ only when the RB 
projection space is the same for all solution subspaces Aii (as above). Interestingly, 
the same additional residual term in the output a also arises when the tensor A{x, •) 
is not symmetric. It is then evaluated with dual cell functions, solutions for a 
problem dual to the cell problems (|2ip llSij . 

3.2 The reduced-basis construction 

Let us choose a sample T> = {xk, 1 < k < p} oi p values for the parameter x in fi. 
Unless some physical properties of the system guides the choice for D, the parameter 
values Xk should be p realizations of a random variable uniformly distributed over 
Q. To accurately solve equation (j2ip for each parameter value x^, we choose an FE 
method with a AA-dimensional FE vector space A^j\f. Typically, M is very large for the 
FE approximations to be accurate. The (n x p) FE approximations [wij^{xk, ■))■ ^ 

for (^uii{xk,-))i k span an (n x p)-dimensional vector space Xnxp C X that contains 
the approximation 

= {ww{xk,y),Xk G V} 

of the solution subspace Ai. 

Remark 3 Computing the (n x p) FE approximations {wij\f{xk, ■))■ ^ can become a 
cumbersome preliminary task to the reduced basis construction when the correspond- 
ing FE matrices are difficult to assemble. That is why the RB method also includes 
some pre-processing to assemble fast those FE matrices. Such a pre-processing is 
very simple in the case where the parametrization of the oscillating coefficients is 
affine (this terminology will be made clearer in section W^) . It might be more diffi- 
cult in other more general cases. In the present work, we only treat the affine case. 
Some more elaborate results in non-affine cases, that are based on the extrapolation 
method introduced in |3 [71) / for instance, will appear in 

First, in the offline stage of the RB approach, we would like to build a A^- 
dimensional RB projection subspace Xn C X that also contains a very close ap- 
proximation of M-^ , thus of M.. Xn will be spanned by a reduced basis {ij)i<j<N 
made of vectors of Xnxp, with N < nxp. Moreover, N <^ M should be sufficiently 
small for the model order reduction to allow a significant gain of computation time. 

Then, in the online stage, for any x £ Q, Wi{x,y) is to be approximated, using 
the RB method, by some vector Wi]\f{x,y) in the Galerkin projection space Xn of 
size that writes 

N 

WiN{x,y) = '^WiNj{x)Cj{y) ■ 

i=i 

The reduced basis (Cj(y))i<j<Ar of Xn is built in order to best control the ap- 
proximation error for outputs through the a posteriori error bounds derived above. 
This is performed in the offline stage as follows. 
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Algorithm 2 (Offline algorithm) We build a reduced basis ('^j(y))i<j<jv from 
Span{wi{xk, ■), Xk £ V,! 1^ i 1^ ^} follows: 

1. for some couple {k^{\),i^{l)) , 1 < k^{l) < p and 1 < i^{l) < n, compute the 
accurate FE approximation tt;jO(x)A/'(3^fcO(i)) y) fof w^iO(i)(xfcO(i), y), element of 

= {wiM{xk,y),Xk &V,l <i <n}, 

2. set 3 = I, ii{y) - 



\wiM{xk,-)\\x' 
3. while j < N, 

(a) compute for every G P and 1 < i < n the (n x p) RB approximations 
Wij{xk,y) G Xj = spanj^fc, 1 < A; < j} for the n cell problems (fl3|). 

(b) for (A;°(j + 1), i°(j + 1)) = argmax ^J^^^J^^^' ^ compute the ac- 

l<k<p,l<i<n \\''^ij[Xk,-)\\x 

curate FE approximation u;iO(i+i)Ar(a;fcO(j+i), y) for WiO(^j^i){xf,o(_j^^,y), 
element of = {wij\f{xk,y),Xk £ V,! < i < n}, 

(c) set ^j+i(y) = - — — where Rj+i is the remainder of the projection 

\\Rj+ny)\\x 

on the j- dimensional reduced basis, 

j 

Rj+i{y) = WiO(^j+i){xkO(^j+i),y) - ^(?i'iO(j-+i)(j;fcO(j+i), ■),Ck)xCk{y), 

k=l 

(d) doj = j + l. 

3.3 Convergence of the reduced-basis method for the cell problem 

The a priori convergence of Galerkin approximations for solutions of continuous 
and coercive elliptic equations like (PT|) is classical. It usually relies on the following 
lemma (see e.g. [19] for a proof). 

Lemma 1 (Cea Lemma) For any 1 < i < n, let Wi be the solution of h21\) and 
WiN its approximation in some N -dimensional Galerkin projection space Xn C X. 
Then we have, for any x € 



\\wi{x,y) -WiN{x,y)\\x < inf \\wi{x,y) - w{y)\\x. 

To conclude that RB approximations like WiN G Xn cl priori converge to vui £ X 
when — > oo, it is then usual to use Lemma [T] in order to a priori prove the 
convergence of the approximation method. 

Lemma 2 // there exists a dense separable subspace V of M and an application 
rjv : V ^ Xn such that 

lim llw - riv(u)|| = 0,Vi; G V , (36) 

then, by Cea Lemma, for any x G and 1 < i < n, RB approximations WiN{x,-) 
converge to Wi{x, •) in the following sense 

lira \\wi{x,y) - WiN{x,y)\\x ■ (37) 

N^oo 
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That is, for all e > 0, there exists a positive integer N{e) such that, Wx £ Q and 
l<i<N, 

\\wiix, •) - WiNix, < e, ViV > N{e) . (38) 

Let us then naturally choose V = ^A, and r^v as the projection operator from 
to Xn for the inner product in X- Unfortunately, the convergence assumed in 
(|36|) can only be shown insofar as we have information about 2?, which amounts to 
knowing how the parameter values are selected to build Xn as N increases. Such 
an assumption is unrealistic since, to choose the right parameter values x^ for P, 
one should already know A4 or some spectral representation of it [16]. So the scope 
of Lemma [2] seems strongly limited, as any a prioi analysis of the RB method in 
general. 

As a matter fact, the RB method is a practical method of order reduction and 
can only be a posteriori shown to converge using reliable and computationally un- 
expensive error bounds that can be evaluated along the RB approximations. 

Note last that, by definition, the Galerkin projection space Xn is built to con- 
verge to the manifold Mj^f = {wi/^{x, y),x G induced by the FE approximations 
Wij\f{x, y) in the sense that, for some given parameter x E Q and 1 < i < ra, there 
exists for all e > a positive integer Mi{e,x) such that 

\\wi{x,y) - Wi^f{x,y)\\x < e,VA/' > ^(e,^) . 

So, let us assume that we have an error estimate for i'^f w{y)£XnxN y) ~ '^iy)\\x 

that is global in parameter space O like in Lemma El in the limit N ^ oo. Even 
then, on account of the pointwise convergence of FE approximations in parameter 
space, the RB method can only converge in the following sense 

lim lim \\wi{x,y) - WiNix,y)\\x = , (39) 

Af^oo A'-*oo 

where Wii\i{x, •) implicitly depends on J\f and where the limits for N and M are not 
reversible. Yet, if the error estimate is also global in parameter space Q with respect 
to the limit — > oo, then the limit for N and M be inverted. 

3.4 Error estimate for the asymptotic homogenized solution 

In the frame of the two-scale homogenization strategy, the asymptotic homoge- 
nized approximation for u'^{x) in the limit e — > is 

uo{x) + e ui (^x, = u*{x) + e Wi (^x, diU*{x) , 

which strongly converges to u'^(x) in (^2) when e — > if G W'^'°°{Q). 

In this approximation, the homogenized solution u* is the solution to the varia- 
tional formulation pU]) of the homogenized equation ^ 

[ A*Vu*-Vv= [ fv+ [ vyv£Hl^{n). (40) 

But in practice, one can only compute an RB approximation for A* , namely 
with entries 

[aUx))^^ = (^I^A{x,y)dyy - sf^ix) , 
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which should be taken into account to estimate the approximation error for the 
asymptotic homogenized approximation. The following lemma ([3]) will show how 
the a posteriori control of the RB approximation allows to control the approximation 
error for the asymptotic homogenized approximation with an RB approach. 

Let us then define an approximation for the asymptotic homogenized ap- 
proximation, 

,{x)+e ^ WiN (x,-) diU*^{x) , 



u 



l<i<n 



where Wijsr is the RB approximation for Wi defined in previous sections, and is an 
approximation for u* that is solution in W^^^^ for the discrete variational problem 



fv 



(41) 



with Whf^^^ C -?/p^(r2) a discrete FE Galerkin projection space associated with a 
mesh of size /i/iom for 0,. We have the following result. 

Lemma 3 Assume that Td is a measurable subset of dQ with a positive [n — 1) — 
dimensional measure (when n > 1) so that a Poincare inequality holds for elements 
of the Sobolev space H^^{i^) = {v £ H^{Q),v |r£,= 0}. 

If the approximations WiN{x, •) converge to Wi{x, •) for all parameter values x in 
n, in the sense 



lim max I esssup\\wii{x,y) - wiiNix,y)\\x > =0 , 



(42) 



then the asymptotic homogenized approximation converges to u*, and so does 
the approximation for the asymptotic homogenized approximation of . That 
is, we have the two results 



and 



lim lim 

N — >oo e^O 



^Im^ lim Ik* (2;) -Un{x)\\l2(q) =0 



u*ix) - u%ix) + e (x, f ) 8.u*ix) - u,^, (x, f ) d.u^,(x, 

l<i<n 



where the two successive limits cannot be inverted. 



Remark 4 As explained in section [373[ the assumption (I42|) can barely be satisfied 
a priori. But in practice, the error bounds derived in the a posteriori analysis of 
section Wl\ allow to check this assumption. The numerical results of Section\^ even 
show that the convergence with respect to N in (j42p is exponential. 



Proof of Lemma\^ To show this result, let us define two quantities, 



E 



N 



u*{x) -u*^{x) 



and 



^Vu (x) ^ y^(^y^-^ _ u%){x) + ^ (VyWi (^x, diu*{x) - VyWiN (^x, diu%{x)^ 

i=l 
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The approximation errors for the asymptotic and homogenized approxi- 
mation of u^{x) now respectively write 



|n*(x) -u^(x)||i2(n) 



E 



u* (x) 
N 



and 



l<i<n 



(x) - 'u^(x) + e ^ \^^Wi (x,-j diU*{x) - WiN (x,-) diU%{x) 



//i(n) 



E 



u*{x) 
N 



L2(Q) 



+ 



E 



N 



+ O (e) 

L2(n) e^o' 



Thus, the proof consists of the two successive results 



lim lim 

N >ooe-+0 



E 



N 



and 



lim lim 

N >ooe^0 



rpSJu* 



L2{Q) 



L2(n) 



. 



(43) 
(44) 



First, let us begin with properties of the homogenized tensor. On account of 
definition (jlSp . A*[x) is a positive definite and continuous matrix. 
Indeed, for any x A*[x) is positive definite 



< a^u • u < aA 


f 

u ■ u + 


I 









^Y^Ui^yWi{x,y)\ 



dy < A*{x)u -u^yuG M" 



since Wi{x, •) is periodic. 

And there exists a positive constant 7*(x) such that 7*(x) is a continuity bound 
for J*(x) 



A*{x)u ■ u < [ u ■ u + 
since the bilinear form in M"" x 

{u,v) 



'^Ui[VyWi{x,y)] 



i=l 



< 'y*{x)u ■ u,yu £ 



|^^ni[Vyii;i(x,y)]j . f ^ i;j[Vj,u;j(x, y)] j dy 



is clearly continuous. 

Moreover, we have uniform continuity because is bounded. That is, there 
exists a real number 7^* > such that, for any x in 57, 7*(x) < 7a*- 

Second, u* G ^hhom ^ ^h^om satisfy variational formulations ([lO]) and 

(jlT]) . We then have the following equality 



A'Vu* ■ Vv 



fv+ V 



that we rewrite with v = (n* — u^) 

[ A'Viu' - u%) ■ V{u' - n^) = / (^^ - A')Vu% ■ V(n* - u'^, 
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Because of the coercivity of A*(x), we finally have the inequality 

qa||V('u* - n^)||i2(Q) < - ^*||oo||Vu^||i2(Q) . 
Moreover, the Poincare inequality for u* — in H^^{Q) writes as follows, 

Wu" - UNh^n) < V\\V{u'' - n^)||L2(Q), 

with a certain constant V which only depends on 0. We have established an error 
estimate for \\Efj 11^2(^2). 

Next, since A''{x) is a positive definite matrix for any x in J7, we deduce the 
following inequality 



Vu* ||2 



X 



< A X,- E 



N 



dx . 



In the limit e ^ 0, on account of the periodicity of A[x, •) and 'VyWiN{x, •), the 
previous inequality rewrites 

2 

L2(f7) 



lim 

e-»0 



rp'Vu* 



< 



J J ^^aA^^ [y^S^i + VyWi{x,y))diU*{x) - {a + VyWiN{x,y))diU*^{x) 



dydx . 



Last, the definition (llSp of the homogenized tensor A* allows to rewrite the 
expression 

A{x,y)^'^{ei + VyWi{x,y))diU*{x) - Vy{ei + WiN{x,y))diU%{x) dy 

i=l 



Y 



and we finally get the following error estimate 



lim 



E 



N 



A* -A 



Iq aA 

the superior bound of which is itself superiorly bounded by 



OA 



+ —\\A* -A 

UA 



* l|2 



Ar||oo||Vn^ll^2(f^) . 



In the end, we have the following error estimates 



lim 



E 



Vu* 

N 



1 



L2(n) a A 

2 



OA 



iVn 



* l|2 



Afllooll V "'7vllL2(Q) 



lim 



E 



N 



L2(o) 



En 



V 



< —\\A* -A 



IVn 



* Il2 



(45) 

i2(o)-a^" 7Vlloo||v«jvllL2(n) • (46) 

They show that the asymptotic homogenized approximations converge if the 
approximate homogenized tensor A^ converges to A*. 

Now, recall that the homogenized tensor A^ converges to A* if the approxima- 
tions Wiiy{x, •) converge to the cell functions Wi{x,-) since we have already obtained 
the following error estimate 



\A* — A^||[^cx)(Q)]nxn 



max {esssup | Sij{x) — sfj{x) |} 
< 77^ max <^ esssup \\wi{x,y) - WiN{x,y)\\x 

l<i<n [ ^izQ 

to derive error bounds for the output s. This concludes the proof of Lemma [3l O 
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3.5 Practical influence of the parametrization 

This section is devoted to the pre-processing used by the RB method in order to 
fast assemble the FE and RB matrices corresponding to the projections of the vari- 
ational formulation (I2ip of the cell problem on the discrete FE and RB Galerkin 
approximation spaces. 

Indeed, for a given family A{x, y) of tensors and a given range Q for parameter 
X values, the solution subspace M for cell problems ()2ip is completely determined 
and fixed. Then, from the theoretical point of view, the way functions w in A4 
explicitly depend on some parameter x G 0, which we call the parametrization, 
should not influence the efficiency of the RB method as a model order reduction 
technique, however it induces the solution subspace Ai. But in practice, the explicit 
parametrization of can significantly account for the efficiency of the RB method, 
because it greatly influences the practical assembling of the matrix and vectors in 
the Galerkin projection method. 

Only piecewise-affine parametrizations (according to the terminology explained 
hereafter) are treated in this work, which allows for a fast, very accurate and simple 
pre-processing of the FE and RB matrices. But the RB method also adapts to 
other types of parametrizations (remember that we already refered to [7] for more 
elaborate results in non-affine cases, based on the extrapolation method introduced 
in ^ jl2] for instance) . 

In the case of an affine parametrization, the assembling of the matrix and vectors 
corresponding to the Galerkin projection of cell problems is always fast and easy. 
By affine parametrization of the cell problems, we mean that A{x, ■) depends on the 
parametrization in an affine manner as follows: for any x S 0, 

_ z _ 

I(x, y) = Ao{y) + e'ix)Ag{y),yy G Y (47) 

9=1 

where: 

- the matrix ^o(y) defines a parameter-independent continuous and coercive bilinear 

form in X X X, 

ao{u,v) = / Ao{y)Vu{y) ■ Vv{y)dy,y{u,v) e X x A! , 

JY 

- the functions G'^ : $7 ^ M are parameter-dependent coefficient functions and 

- the matrices Aq{y) define parameter-independent continuous bilinear forms in 

aq{u,v) = J^Ag{y)Vu{y) ■Vv{y)dyy{u,v) €XxX ■ 

With such affine parametrizations, the numerical RB treatment of cell problems 
is straightforward. Let us detail its implementation. 

First, we follow the offline algorithm presented in the section 13.21 At each step 
of the offline stage, a cell problem ()14p for some parameter value in P = {xj., 1 < 
k < p} is to be explicitly solved in order to build the reduced basis (Ci(y))i<j<jv 
span{tt)(x, y), X G V}- We choose to use an FE method for this initial step of the 
offline stage that consists of accurately solving the variational formulation (|2ip with 
conforming Pi Lagrange finite elements and a fine mesh for Y. 
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The solution space X is discretized into the vector space Xj\f of continuous, 
piecewise hnear functions. Let 7y be a conformal mesh for the n-torus Y = [0, 1]" 
made of A/j elements {'^k)i<k<Aft of size hy- We write (pk the FE basis functions 
associated with the M nodes yk in Y. Now, for < q < Z,we define the FE matrices 



Mq e M-^ with entries 



/I] 

for any 1 < «,j < AA, and n FE data vectors F^^/ G M-^ (1 < ^ < with entries, for 
any 1 < i < Af, 

{Fk,i)i = / M{y)ei ■ Vy(t)i{y)dy . 



JY 

Then, for any x in P, 1 < i < n, we compute the FE approximate solution 

M 

ww{x,y) E Xjv = '^WiXk{x)4'k{y) 
k=l 

for the cell problem (j21|) . that satisfies 

Mo + Yl ^'i^Wq ) wMx, yi) = (pQ^i + (48) 

9=1 / \ fc=l J I 

at each node y;, 1 < / < N . The FE problem (|48|) can then be fast and very 
accurately assembled through linear combinations of the matrices (-Mg)o<g<z and 
vectors (i^g, 2)o<A:<Zi 1 < ^ < 'i-, that are to be kept into memory. 

Then, in the online stage, we would like to treat fast cell problems (I2ip for any 
parameter value x G O. Let us project the FE matrices and vectors on the RB space 

^RB ^ ^^Mg^, 0<q<Z 

and 

with ^ the A^x (nN) matrix with columns ^fc, 1 < A; < n x A^. The RB approximation 
WiNix,y) = J2k=i'^iNk{x)Ckiy) for cell function w{x,y) is solution of the linear 
system 

M(f ^ + E ©'(^X"" I yi) = f ^oT + E 0'(^)^9T1 ' (49) 

9=1 / V fc=l /; 

that is fast assembled through linear combinations in the present affine case. And, 
for 1 < i,j < n, the outputs are easily given by 

nN / Z \ 

1=1 \ q=l J ^ 

Moreover, an error bound for each RB approximation can also be derived fast in the 
online stage, following an offline-online strategy similar to that applied to the RB 
output prediction \TE\. 

So, affine parametrization obviously allows for a fast assembling of the RB matrix 
and vectors in (j49|) . It can be considered as an ideal frame for an efficient RB method, 
because the CPU time for the online solution of one cell problem (I49p actually scales 
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like the CPU time for solving a linear system of size A^, and because the offline stage 
is actually very short in comparison with the online computation of a large number of 
cell functions. The possibility of a similar gain of computation time is not as obvious 
in the case of non-affine parametrizations, when (|47p is not valid anymore. Then, 
the assembling of matrices and vectors needs more elaborate techniques [31IT2]. Yet, 
an efficient RB approach is still possible for quite a few situations, as shown in [7|. 

We finally treat the case of piecewise affine parametrizations that can be re- 
casted into the class of affine parametrizations. The pre-processing that we propose 
in this case lies on the fact that, in practice, the RB approximations are numerical 
approximations for FE approximations, and not for the "true" cell functions. Then, 
to apply the RB method to the parametrized FE approximations, it is possible to 
consider a "global" parametrization made of a parameter for the oscillating coeffi- 
cients and of another parameter for the FE method (for instance, the geometrical 
features of the mesh). 

We deal with oscillating coefficients A{x,-) parametrized in a piecewise affine 
manner as follows: 

- for each x £ Q, the cell Y can be partitioned into d non-overlapping Yk{x) open 

subsets (d S N* should be fixed) such that Y C Ufc=i '^k{x), 

- there exists d non-overlapping reference open subsets Y^ such that Y C ULi ^fc ' 

- for each x £ Q, there exists d affine homeomorphisms, <I>fc(x,-) : Y^ Yk{x), 

l<k<d, 

- and for every 1 < k < d, the family of functions ( A{x,^{x, ■))] restricted to 

Y^ can be parametrized in an affine manner as defined in (j47p by 

_ z _ 

I(x,$(x,y)) =io(y) + 5]e'?(x)i,(y),Vy G Y^ . (50) 

q=l 



The function <I>(x, •), defined almost everywhere in Y by 

^{x,y) = ^k{x,y),yy eY^,l<k< d, 

maps a "reference" cell onto the cell with parameter value x. After the mapping, 
the family of cell problems defined with these piecewise affine oscillating coefficients 
can then be treated as if the parametrization was affine like in ()47p . provided one 
take into account the stretching of the domain at each parameter value x. 
For this, we define 2(Z + 1) tensors of rank 3, ( I and ( 

V /l<fe<Z+l V / l<fe<Z+l 

by: 

Aft N M 



^^ = Y^Y^Y^\\ My)'^yMy) ■ '^y<Pj{y)dy ) ® ® ej, (51) 
1=1 i=i j=i ^-''^i ' 

Aft N N ^ ^ _ X 

^k = Yl^^( Ak{y)ei-Vy(^j{y)dy] ei(S)ei(^ej . (52) 



1=1 i=l j=l 



An accurate pre-processing in the piecewise affine cases is then possible that assem- 
bles fast Pi-FE matrices, and corresponding RB matrices, by adding a mapping step 
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to the linear combinations of the affine cases. Namely, with the family of vectors 
that accounts for the stretching of the mesh elements, 

M 

V(x) = det ((Vj,</.(x, ) , Vx e 0, 
1=1 

we easily get the FE matrix for any parameter value x in through the formula 

z 



V- Mo + ^e'?(x)M< 



9=1 



and so on for the RB matrix. Note also that the reduced basis should then be 
orthonormalized at each parameter value x in fi, because the inner product matrix 
also changes for each x. 



4 Numerical results 

We now show numerical results for the RB approximation of a seemingly non-affine 
two-dimensional problem that is brought back to the affine setting after mapping 
of the cell Y. We do not show the MP-RB treatment of more general piecewise 
continuous parametrizations in this work, but elementary results for one dimen- 
sional problems can be found in [Tj. The two-dimensional problem is chosen here to 
show the efficiency of the RB method in a classical situation for the homogenization 
theory. To fix ideas, it consists of homogenizing the conductivity of a heteroge- 
neous composite material in a domain Q, where a two-dimensional matrix is full of 
inclusions with varying positions and conductivity properties. 



4.1 Definition of the problem 



2/2 

1 

0.75- 



0.25 










Qo 











0.25 



y2 
1 













Q{x) 











0.75 1 yi 



h 



ci 



1 yi 



Figure 1: For each parameter value x, the cell with inclusion Q{x) (on the right) is 
mapped through the piecewise affine homoemorphism $(x, •) from a reference cell 
with inclusion Qo (on the left). 

For n = 2 and / = 0, we supply the problem ([1]) with the mixed boundary 
conditions 



(BC) 



U'{1,X2) 





+1 



u^{xi, 1) 
^^Vn^ • n\ 



(xi.O) 



(53) 



We define at each point x G $7 a single rectangular inclusion Q{x) C y in the 
ceh Y = [0,1]2, Q(x) = {{yi,y2)\0 < h{x) < yi < q(x) < 1,1 < i < 2} (FigH]). We 
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also write I2 the second-order identity tensor and 1q(x) the Q{x)-test function, such 
that, for every y £Y, lQ(x)(y) is one if y G Q{x) and zero otherwise. 

For ah X G 17, the oscillating coefficients A'^ix) = I2 + ^i(x,e~^x) are locally 
periodic with a Y-periodic function Ai{x, y) = 6'(x)1q(x)(2/)-^2 that is constant inside 
and outside of the inclusion Q{x). 

We want to homogenize the problem ([T]) associated with oscillating coefficients 
parametrized by the multiparameter (61, ci, 627 C2, ^)(a;), that is function of x G i7 
and takes value in [.25-5; .25 + 5]2 x [.75-5; .75 + (5]2 x [-6^°; 0], where J g]0; .25[ and 
6^ G]0; 1[. For the FE matrices to be easily assembled, we define a "reference" cell 
problem with a centered inclusion Qq = [0.25; 0.75]^ (Fig{T]). Then, at each point 
x G fi, the inclusion Q{x) can be mapped on Qq as explained in section [331 



4.2 Offline computations 




Figure 2: Maximal relative error bound max ^ — ^ — for the RB approx- 

l<i<2,Xk&V \\Wi{Xk,-)\\x 

imations WiN^Xk, •) of the initial sample used for the RB construction (left picture), 
and for the RB approximations WiN{zk, •) of a test sample G A (right picture), in 
log-scale with respect to the size of the growing reduced basis. 

A FE approach is developped for mapped cell problems in Y with the "reference" 
inclusion Qq. More precisely, we use classical Pi simplicial Lagrange finite elements 
on a quadrangular, uniform and affine FE mesh, divided in isosceles triangles with 
base along direction y2 = —yi and size hy in each direction ei and 62- The mesh 
is fixed and adapted to the "reference" domain in the sense that the boundaries of 
the inclusion Qq are multiples of hy- 

We choose a random initial sample T? of parameter values that is uniformly 
distributed over the multiparameter range. A reduced basis is then built for any 
parameter point x after mapping with <I>(x, •) the solutions Wi^Xk, ^{xk,y)) selected 
by the offline algorithm of section 13.21 for an initial sample P of p = 50 parameter 
values. Numerical results are shown forh 6 = .1, 9^ = .99 and hy = .1 in Figures 
[2] and [3l Note that the contrast between the coefficients inside and outside the 
inclusions can grow up to 1/100, which is makes our expremient quite a stringent 
test. 

The relative a posteriori error bounds for the RB approximations at the param- 
eter values of the initial sample are computed at each step of the offline algorithm. 
The maximal error bound in this initial sample decreases exponentially with the 
size N of the reduced basis (Fig. [2]). The effectivity of the a posteriori estimation 
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^<i,j<2,ZkeA \Sij{Zk)\ 

a random test sample A of parameter values in 0. 



is checked all along the RB construction (we found r]j^{wi{xk, •)) G]1.4;3.5[ for all 
1 < k < p and 1 < i < n in the numerical experiment corresponding to Fig. [5]). 
Note that the offline algorithm selects (almost always alternatively) cell functions 
for the both cell problems, in direction ei and 62- Then, for N = 2, one cell function 
per direction only spans the reduced-basis, which strongly amplifies the RB approxi- 
mation errors for the cell problem corresponding to the second direction represented 
in the reduced basis. 

The reduced basis is then tested for another sample A = < A; < p} of 

parameter values in Q. The maximal a posteriori error in this test sample still 
decreases very fast (exponentially with the size A'' of the reduced basis), but the 
rate of decrease is slightly smaller than that of the initial sample used for the RB 
construction (Fig. [2]). This shows that the initial sample P was not an optimal 
choice to compute a reduced basis for any x £ yet it still allows for efficient 
Galerkin approximations with any A. 

Besides, the actual output approximation error for this test sample scales as the 
square of the actual approximation error for cell functions (see Fig. [Sj obtained in 
the same numerical experiment than Fig. That is, the RB approximations are all 
the more efficient for the outputs, and the approximation errors scale like the error 
bounds derived in section [3TT1 The effectivities of the error bounds of section [3TT] are 
indeed hardly bigger than one (we found r]N{wi{zk, •)) G]1.3;3.9[, for all 1 < k < p 
and 1 < i < n, in the numerical experiment corresponding to Fig. [3|0. 

4.3 Online computations 

After building a reduced basis with the greedy algorithm from the previous FE 
approximations, we use the RB method to compute online RB approximations for 
cell functions as linear combinations of the RB basis functions. For this online 
stage, we develop an FE method for the homogenized problem ^ and use classical 
Pi simplicial Lagrange finite elements on a quadrangular, uniform and affine FE 

*Note that the maximal relative a posteriori error bound and the maximal actual error in Figures 
[5] and [3] are not obtained for the same parameter value Zk, hence the discrepancy between their 
ratio and the effectivity riN{wi{zk, ■)) we measured. 
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ratio 

N 1 Kf 
l\ 1 J\ 


VP — ouj 

IIY 


digUilLlllli 


f ^hom — 3 ^ 

6 


ix-D lor -^j^ 


E Hi lOr 


1/0 


1 


1 7 


2E-'^ 




97 ^ 


1/5 


1^-1 


15 s 


2E-^ 


410+330 = 740 s 


3100 s 


1/20 




42 s 


2E--' 


16+10 = 26 s 


520 s 


1/20 




53 s 


2E-'-^ 


1600+1000 = 2600 s 


37000 s 



Table 1: CPU time (in seconds) needed by a Matlab code with an Intel Pentium 
IV processor (3.0 GHz/1 Go) to approximate the FE matrix for the homogenized 
problem either with a direct FE approach or with an RB method. In the RB ap- 
proach, one has to take into account the RB construction (offline algorithm with 
a sample of p parameter values), the online computation of one homogenized so- 
lution, plus possibly the online a posteriori estimation error (hence the two terms, 
solution+estimation, in the RB online column). 



( ^hom — 3 A 

V e ~ 2) 


\\u'^ — U*\\l2 


W'^Af ~ ""iV HI 


l|Vr,||i2 


\\^y{WiJ\f - WiN)\\L2 


(theory) 


< Cie 




<C2V^ 




hy = lE-^ 


(e = 2.0^-^) 


1.2^-* 


{^e = lAE-') 


2.9E-'^ 


hy = lE-^ 


(e = 2.0^-3) 


4.7^-3 


iV~e = 4.5E-'') 


l.OE-'^ 


hy = 5E-'^ 


(e = 2.0^-^) 


S.IE-'-^ 


{^=1AE-') 


8.6E-^ 


hy = 5S-2 


(e = 2.0E-'^) 


I.IE-'-^ 


iV~e = 4.5E-^) 


3.0^-2 



Table 2: Theoretical correction error for the homogenized solution, and RB numeri- 
cal approximation error for the homogenized solution when 6 = .2, 6^ = .99, p = 50 
and N = 20. 



mesh divided in isosceles triangles with base along direction X2 = —xi and size hhom 
in each direction ei and 62- 

The RB computations are performed in the step of the numerical homogenization 
strategy where the values of the homogenized coefficients are collected, as outputs 
of the cell functions, at some quadrature points in O that are necessary for the 
computations of the entries of the FE matrix in the homogenized problem ([8]). The 
CPU time needed for computing these outputs is compared between the RB and FE 
methods (Tab. [T]), where the RB method includes an online a posteriori estimation 
of its approximate solution. 

In calculations of Table [H the RB method has been applied with a reduced basis 
of size = 20 starting from an initial parameter sample D of size p = 50. The 
main result of Table [1] is that the ratio of the RB computation time on the FE 
computation time scales like N/M, the ratio of the numbers of degrees of liberty in 
the RB and (direct) FE methods. 

So we can distinguish between two main regimes. The more unfavourable regime 
is the case of large ratios N/M, which corresponds to cases where one need only small 
precision for the correction (jlip (large hy). Then, the RB method is likely to be 
faster than a direct FE method in the frame of many queries of the homogenized 
solutions. Note that in such a situation, the computation time spent by the offline 
algorithm is not even an issue. It is then possible to enlarge the initial parameter 
sample 2? (take a larger p). This increases the computation time spent by the offline 
algorithm in the RB construction but improves the quality of the RB approxima- 
tions. 

On the contrary, the favourable regime corresponds to small ratios N/J\f, where 
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the correction is sought very accurate (smaU /ly). Then, the numerical results for 
the RB approximations of the cell functions show that there is an important gain of 
computation time, while there is no significant loss of numerical precision (Tab. [2]). 

5 Conclusion and perspectives 

We have shown in the present work that, for a prototypical class of parametrized cell 
problems (with piecewise affine oscillating coefficients), the reduced-basis approach 
applies and significantly reduces the time needed to compute a large number of 
parametrized cell problems in homogenization, in comparison with an FE method. 

Some interesting questions concerning the extension of the RB approach in ho- 
mogenization remain, mainly linked to the treatment of a larger class of parametrized 
cell problems: in particular, other geometries for more realistic cell problems should 
now be addressed, other boundary conditions for the cell problems (including the 
treatment of oversampling techniques), less regular oscillating coefficients (with 
many inclusions in varying amount). Also, the same questions as those examined 
in the present work for scalar elliptic equations could be asked for the Stokes-Darcy 
equations in porous media, or for the equations of linear elasticity in two- and 
three-dimensional contexts. Further developments of the RB methodology are then 
needed that may lead to interesting (fast) approaches in homogenization. Since a 
major issue in homogenization is the limitation of the time computation, speed- 
ing up the homogenization procedures could inevitably bring new possibilities of 
refinements (perhaps like reiterated oversamplings to improve the accuracy of the 
corrector term). 

In any case, we believe that our result is interesting in the frame of many of 
the commonly used homogenization strategies, namely all those that ask for solving 
a computationally demanding number of parametrized cell problems. This is true 
provided the type of the parametrization can be handled with our RB approach. 
Among those homogenization strategies, the two-scale homogenization strategy is 
well known and much used in practice. That is why we have chosen this frame for 
our numerical experiments. But other homogenization strategies, which are used for 
non-locally-periodic oscillating coefficients, can also be treated with an RB approach. 

For example, stochastic homogenization also asks for solving a large number 
of parametrized cell problems in the frame of local approximations of the homog- 
enized tensor [B]. The homogenization of locally deformed oscillating coefficients, 
A''{x) = A{^~^{e~^x)) with ^ a diffeomorphism, in the frame of deterministic ho- 
mogenization [8] or of stochastic homogenization [5], A''(x) = A{^~^ {€~^ x , uj)) with 
Lo an element of a probability space, by nature, also demand for solving parametrized 
cell problems. 

More general cases, often computationally demanding, also rely on the compu- 
tation of a large number of cell problems, and offer a frame for an application of 
the RB approach. Among those homogenization strategies, the heterogeneous mul- 
tiscale method (HMM), that averages over a large number of cell problems, could 
directly make use of our RB approach when cell problems are correctly parametrized. 
Another one, the multiscale finite-element method (MsFEM), also averages over nu- 
merous cell problems. Yet, the range of geometries for those cell problems is often 
larger, and it is still not obvious that model order reduction techniques may speed 
up the MsFEM computations. 

Although we have not tested all the above mentioned possible improvement, we 
believe that our work is likely to improve a large number of existing homogenization 
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strategies. Definite conclusions on the validity of our approach in such settings will 
hopefully be obtained soon. 
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